MODULE transport

    IMPLICIT NONE


CONTAINS

!------------------------------------------------------------
!Transport routine to call subroutines
!------------------------------------------------------------
SUBROUTINE aetransport(pos)
  
 USE definitions
 USE deposition
 USE erosion
 USE avalanche_modules
 USE functions
 
 IMPLICIT NONE
 
 INTEGER, INTENT(in), DIMENSION(2)	:: pos
 INTEGER, DIMENSION(2)			:: npos, bpos
 INTEGER					:: e, d, trigx, trigy
 REAL                                    :: L, erhight
        
 npos = pos 
 
 d = 0
 stackcount = 0
 erhight = erode(npos, 1)
 erbuf = erbuf + erhight
 IF(erhight>0) THEN                                   !Erode slab from cell

    IF(full == 1) THEN
       local_eroded_hight(pos(1),pos(2)) = local_eroded_hight(pos(1),pos(2)) + erhight
    ENDIF

    CALL avalanche(npos)                              !Call avalanche routine
    CALL calcl(npos, L)                               !Generate path length L
    DO WHILE(erhight>0)
       d=d+1
       bpos = npos
       CALL find_next_cell(npos, trigx, trigy, L, erhight)    !Move slab to next cell
       IF(trigx == 1) THEN
          xsinkbuf = xsinkbuf + erhight
          sourcebuf = sourcebuf + erhight
          IF(xsink == 1) THEN
             erhight = 0
             EXIT
          ENDIF
       ENDIF
       IF(trigy == 1) THEN
          sourcebuf = sourcebuf + erhight
          IF(ysink == 1) THEN
             erhight = 0
             EXIT
          ENDIF
       ENDIF
       IF(depose(npos,erhight,L)==1) THEN
          EXIT
       ENDIF
    ENDDO
    trigx = 0
    trigy = 0
    CALL avalanche(npos)
    path(d) = path(d) + 1
    pn = pn + 1
 ENDIF

END SUBROUTINE aetransport
!------------------------------------------------------------

SUBROUTINE calcl(pos, L)

  USE definitions

  IMPLICIT NONE
 
  REAL, INTENT(out) :: L
  INTEGER, INTENT(in), DIMENSION(2) :: pos
  REAL, PARAMETER :: deltaL = 0.2
  REAL :: rand

  CALL RANDOM_NUMBER(rand)
  
  IF(gradL == 0) THEN
     L = L0 
  ELSEIF(gradL == 1) THEN
     L = L0 + (L1-L0)/REAL(ysize)*pos(2)
  ENDIF

END SUBROUTINE calcl

!------------------------------------------------------------
!Routine to find next cell according to wind direction
!----------------------------------------------------------
SUBROUTINE find_next_cell(pos, trigx, trigy, L, mass)

	USE definitions
	USE functions

	IMPLICIT NONE

	INTEGER	:: i, j, k
	INTEGER, INTENT(inout), DIMENSION(2) :: pos
        REAL, INTENT(inout) :: L
        REAL, INTENT(in) :: mass
        INTEGER, DIMENSION(2) :: bpos, rpos
        INTEGER, INTENT(out) :: trigx, trigy
	REAL	::	phi, sump, rand
	REAL, DIMENSION(8)	:: p	

!Set probabilities for propagation to next cell
!Indexing of probabilities according to the following scheme:

 ! 		x-1	x	x+1	
 !  					
 ! y+1	7	4	8	
 ! 					
 ! y	1	C	2	
 !					
 ! y-1	5	3	6	
 !					
        bpos = pos
	phi = ATAN(windy/windx)

	DO i=1,4
		p(i) = windx*relpos(1,i)+windy*relpos(2,i)
		IF(p(i)<0) p(i) = 0
	ENDDO

	DO i=5,8
		p(i) = (windx*relpos(1,i)+windy*relpos(2,i))/2.0
		IF(p(i)<0) p(i) = 0
	ENDDO

!Set sum of probabilities to one

	sump = 0
	DO i=1,8
		sump = sump + p(i)
	ENDDO
	p = p/sump

!find next cell
	CALL RANDOM_NUMBER(rand)
	sump = 0
	DO j=1,8
		IF (rand > sump .AND. rand <= sump + p(j)) THEN
			pos = pos + relpos(:,j)	
		ENDIF
		sump = sump + p(j)
	ENDDO

	IF(pos(2)>ysize .AND. ylock > 0) pos(2) = pos(2)-1
	IF(pos(2)<1 .AND. ylock > 0) pos(2) = pos(2)+1
!Calculate flow on adjancent cells for full diag
        IF(full == 1) THEN
           k = ABS(pos(1)-bpos(1))+ABS(pos(2)-bpos(2))
              rpos = pos
              CALL make_periodic(rpos)
           IF(k == 1) THEN
              local_aeolian_flux(bpos(1),bpos(2)) = local_aeolian_flux(bpos(1),bpos(2)) + 1.0/2.0*mass
              local_aeolian_flux(rpos(1),rpos(2)) = local_aeolian_flux(rpos(1),rpos(2)) + 1.0/2.0*mass
           ELSEIF(k == 2) THEN
              local_aeolian_flux(bpos(1),bpos(2)) = local_aeolian_flux(bpos(1),bpos(2)) + 1.0/2.0*mass
              local_aeolian_flux(rpos(1),rpos(2)) = local_aeolian_flux(rpos(1),rpos(2)) + 1.0/2.0*mass

              local_aeolian_flux(rpos(1),bpos(2)) = local_aeolian_flux(rpos(1),bpos(2)) + 1.0/4.0*mass
              local_aeolian_flux(bpos(1),rpos(2)) = local_aeolian_flux(bpos(1),rpos(2)) + 1.0/4.0*mass
           ENDIF
        ENDIF
!Calculate transport length for full diag
        IF(full == 1) THEN
        local_transport_length(rpos(1),rpos(2)) = local_transport_length(rpos(1),rpos(2)) + &
             (rpos(1)*windx + rpos(2)*windy)/(SQRT(windx*windx + windy*windy))*(1+par2*ATAN(gradient_h(rpos)))
        ENDIF
!Reduce pathlength
        rpos = pos - bpos
        bpos = pos
     	CALL make_periodic(pos)
        L = L - slabsize*(rpos(1)*windx + rpos(2)*windy)/&
             (SQRT(windx*windx + windy*windy))*(1-par2*ATAN(gradient_h(pos)))
!Determine if border is crossed
 
        IF(ABS(bpos(1)-pos(1)) > 0) THEN
           trigx = 1
        ELSEIF(ABS(bpos(1)-pos(1)) == 0) THEN
           trigx = 0
        ENDIF

         IF(ABS(bpos(2)-pos(2)) > 0) THEN
           trigy = 1
        ELSEIF(ABS(bpos(2)-pos(2)) == 0) THEN
           trigy = 0
        ENDIF

END SUBROUTINE find_next_cell
!------------------------------------------------------------

END MODULE transport
